<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta name="generator" content="MATLAB R2016a"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><title>Fines migration, hydrodynamic release, and deposition (formation damage) in a five-spot well pattern - case 1</title><style type="text/css">
* {margin: 0; padding: 0;}
body {text-align: start; line-height: 17.2339992523193px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Consolas, Inconsolata, Menlo, monospace; font-style: normal; font-size: 14px; font-weight: normal; text-decoration: none; white-space: normal; }
h1, h2 {font-weight: normal;}
.content { padding: 30px; }

.S0 { margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S1 { line-height: 26.3999996185303px; min-height: 24px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-size: 22px; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 15px; margin-right: 10px;  }
.S2 { min-height: 0px; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S3 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S4 { min-height: 0px; font-weight: bold; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S5 { line-height: 20.576000213623px; min-height: 20px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-size: 16px; font-weight: bold; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 9px; margin-right: 10px;  }
.S6 { min-height: 0px; font-family: monospace; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S7 { margin-left: 3px; margin-top: 10px; margin-bottom: 10px; margin-right: 3px;  }
.S8 { line-height: 15.5926666259766px; min-height: 18px; white-space: nowrap; font-size: 12.6666669845581px; white-space: nowrap; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S9 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(34, 139, 34); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S10 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S11 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S12 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 10px; margin-bottom: 9px; margin-right: 10px;  }
.S13 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(160, 32, 240); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S14 { font-family: Helvetica, Arial, sans-serif; margin-left: 0px; margin-top: 10px; margin-bottom: 20px; margin-right: 0px;  }
.S15 { text-align: left; line-height: 21px; white-space: pre-wrap; white-space: pre-wrap; margin-left: 56px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S16 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S17 { margin-left: 3px; margin-top: 10px; margin-bottom: 4px; margin-right: 3px;  }
.S18 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); font-size: 12.6666679382324px; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S19 { color: rgb(64, 64, 64); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S20 { min-height: 0px; color: rgb(0, 95, 206); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }

.LineNodeBlock {margin: 10px 0 10px 0;}
.LineNodeBlock+.paragraphNode {margin-top: 10px;}
.lineNode {padding-left: 10px; background-color: #F7F7F7; border-left: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.inlineWrapper:first-child .lineNode,.inlineWrapper.outputs+.inlineWrapper .lineNode {padding-top: 5px; border-top: 1px solid #E9E9E9;}
.inlineWrapper:last-child .lineNode,.inlineWrapper.outputs .lineNode {padding-bottom: 5px; border-bottom: 1px solid #E9E9E9;}
.lineNode .textBox {white-space: pre;}
.outputGroup {margin: 2px 0 2px 0; padding: 2px 2px 2px 4px;}
.outputRegion {}
.outputParagraph {color: rgba(64, 64, 64, 1); padding: 10px 0 6px 17px; background: white; overflow-x: hidden;}
.inlineWrapper:last-child .outputParagraph {border-bottom-left-radius: 4px; border-bottom-right-radius: 4px;}
.outputParagraph:empty {margin: 0;}
.inlineElement .symbolicElement {}
.embeddedOutputsSymbolicElement .MathEquation {margin-top: 5px; margin-bottom: 5px;}
.embeddedOutputsSymbolicElement .MathEquation.displaySymbolicElement {margin-left: 15px;}
.embeddedOutputsSymbolicElement .MathEquation.inlineSymbolicElement {margin-left: 5px;}
.embeddedOutputsSymbolicElement {overflow-x: auto; overflow-y: hidden;}
.embeddedOutputsSymbolicElement { overflow: initial !important;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre; word-wrap: initial; min-height: 18px; max-height: 250px; overflow: auto;}
.textElement {padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement { max-height: none !important; overflow: initial !important;}
.veSpecifier {}
.veContainer {}
.veSpecifierBox {height: 400px; width: 500px;}
.veSpecifier .veTable {padding-top: 3px; padding-bottom: 4px;}
.veSpecifierBox .veSpecifier .veContainer {position: relative; width: 100%; height: 370px;}
.veSpecifierBox .dijitDialogPaneContent {width: 97% !important; height: 88% !important;}
.veSpecifier .veTable .rowHeadersWrapper {padding-bottom: 0;}
.veSpecifier .veTable .scroller .variableEditorRenderers {padding-right: 3px; -webkit-user-select: none; -moz-user-select: none; -ms-user-select: none;}
.veVariableValueSummary {display: inline-block;}
.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {visibility: hidden; z-index: 0;}
.veMetaSummary {font-style: italic;}
.veSpecifier .veTable .scroller {overflow: hidden;}
.veSpecifier .veTable:hover .scroller {overflow: auto;}
.veSpecifier .veVariableName,.veSpecifier .veDimensionFont {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.veSpecifier .veVariableName {padding-top: 3px;}
.veSpecifier .veDimensionFont {font-style: italic; color: #9A9A9A;}
.veSpecifier .scroller::-webkit-scrollbar-track {background-color: white;}
.veSpecifier .scroller::-webkit-scrollbar-corner {background-color: white;}
.veSpecifier .veTable .topRowHeaderWrapper {border: none; background-color: #F8F9FA;}
.veSpecifier .mw_type_ListBase.showCellBorders,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper,.veSpecifier .veTable .verticalScrollSpacer,.veSpecifier .veTable .horizontalScrollSpacer {border: none;}
.veSpecifier .veTable .dataScrollerNode {border: 1px solid #BFBFBF;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode,.veSpecifier .veTable .dataBody {font-family: Arial; font-size: 13px;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode {color: #7F7F7F;}
.veSpecifier .veTable .dataBody {color: #000000;}
.veSpecifier .veTable .columnHeaderNode .cell .drag,.veSpecifier .veTable .columnHeaderNode .cell .header,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {background-color: #F8F9FA;}
.veSpecifier .veTable .columnHeaderNode .cell .dragBorder {border-right: 1px solid #F8F9FA;}
.veSpecifier .veTable .rowHeaderNode .cell {padding-top: 3px; text-align: center; border-bottom: 1px solid #F8F9FA;}
.veSpecifier .veTable .dataBody .cell .plainText {text-align: right;}
.veSpecifier .veTable .dataBody .row .cell {border-bottom: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.embeddedOutputsVariableElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre-wrap; word-wrap: break-word; min-height: 18px; max-height: 250px; overflow: auto;}
.variableElement {}
.embeddedOutputsVariableElement.inlineElement {}
.inlineElement .variableElement {}
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement {min-height: 18px; box-sizing: border-box; font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.matrixElement .variableValue {white-space: pre; display: inline-block; vertical-align: top; overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {padding-top: 0 !important; max-height: 100px;}
.embeddedOutputsMatrixElement.inlineElement .veVariableName {padding-top: 0;}
.inlineElement .matrixElement {max-height: 300px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer {white-space: nowrap; padding-top: 10px; margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide {display: none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis {margin-bottom: -3px;}
.matrixElement { max-height: none !important; }
.dijitTooltipDialog .dijitTooltipContainer .dijitTooltipContents .alertPlugin-alertMessage {min-width: 12px; max-width: 400px; max-height: 300px; overflow: auto;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar {width: 11px; height: 11px;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-track {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-corner {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {border-radius: 7px; background-color: rgba(0, 0, 0, 0.1); border: 2px solid rgba(0, 0, 0, 0); background-clip: padding-box;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0.1);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.alertPlugin-alertLine {position: absolute; display: initial; width: 40px; text-align: right; cursor: text;}
.alertPlugin-onTextLine {visibility: hidden;}
.alertPlugin-hasTooltip .alertPlugin-warningImg,.alertPlugin-hasTooltip .alertPlugin-errorImg {cursor: pointer;}
.alertPlugin-isStale {-webkit-filter: opacity(0.4) grayscale(80%); filter: opacity(0.4) grayscale(80%);}
.alertPlugin-alertLine .alertPlugin-errorElement {display: inline-block; margin-right: 4px;}
.alertPlugin-errorImg {position: relative; display: inline-block; background-color: rgb(204, 55, 41); width: 14px; height: 14px; margin-top: 1px; -webkit-border-radius: 7px;}
.alertPlugin-alertLine .alertPlugin-warningElement {display: inline-block; margin-right: 3px;}
.alertPlugin-warningImg {position: relative; display: inline-block; width: 16px; height: 15px; overflow: hidden;}
.diagnosticMessage-wrapper {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {color: rgb(255,100,0); text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {color: rgb(230,0,0); text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart {white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {white-space: pre;}
.embeddedOutputsWarningElement{min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsErrorElement {min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.alertPlugin-warningImg {background-image: url()}
.alertPlugin-errorImg {background-image: url()}
</style></head><body><div class = "content"><div class = 'SectionBlock containment active'><h1 class = "S1"><span class = "S2">Fines migration, hydrodynamic release, and deposition (formation damage) in a five-spot well pattern - case 1</span></h1><p class = "S3"><span class = "S2">In this demonstration example we will reproduce the two-dimensional problem for formation damage caused by massive injection of a CO2-particle mixture as discussed by </span><span class = "S4">Sbai and Azaroual (2011) </span><span class = "S2">in their example described in section 5.1. This involves pumpage of a salt aqueous brine from a 3000m deep reservoir through four square corner's wells with an equal injection rate. CO2 phase is injected from a well at the center of the domain with a total injection flow rate of 1.15 x 10^6 tons/year. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S5"><span class = "S2">Reservoir, fluids,and particulate suspensions properties setup </span></h2><p class = "S3"><span class = "S2">Let's first specify the simple reservoir geometry which is a one layer of 500m x 500m x 5m along each space dimension. Additionally, we will explicitly store some required variables in the '</span><span class = "S6">Grid</span><span class = "S2">' structure as reequired for next computational tasks. These are: </span><span class = "S6">Nx</span><span class = "S2">, </span><span class = "S6">Ny</span><span class = "S2">, </span><span class = "S6">Nz</span><span class = "S2"> which corresponds to the number of cells along each space dismension, respectively. N: the total number of grid cells, </span><span class = "S6">hx</span><span class = "S2">, </span><span class = "S6">hy</span><span class = "S2">, </span><span class = "S6">hz</span><span class = "S2">, which are the uniform spacing along each space dimension, </span><span class = "S6">V</span><span class = "S2">: the volume of grid cells, </span><span class = "S6">K</span><span class = "S2">: the initial distribution of the permeability, </span><span class = "S6">compr</span><span class = "S2"> the total compressibility factor, </span><span class = "S6">por</span><span class = "S2">: the initial distribution of the porosity. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Domain size along X, Y &amp; Z directions</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">dom = RecDomain([500,500,5]);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Dx = dom.Lx; Dy = dom.Ly; Dz = dom.Lz;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.Nx = 101;  Grid.Ny = 101;   Grid.Nz = 1; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Nx = Grid.Nx;   Ny = Grid.Ny;    Nz = Grid.Nz; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hx = (Dx/Nx);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hy = (Dy/Ny);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.hz = (Dz/Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.N = Grid.Nx*Grid.Ny*Grid.Nz;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">N = Grid.N;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);                        </span><span class = "S9">% grid cells volumes</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.K = 0.85e-12.*ones(3,Nx,Ny,Nz);                      </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">Grid.compr = 4.4e-10.*ones(Grid.Nx,Grid.Ny,Grid.Nz);     </span><span class = "S9">% compressibility</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">Grid.por   = 0.3.*ones(Grid.Nx,Grid.Ny,Grid.Nz);         </span><span class = "S9">% porosity</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div></div><p class = "S12"><span class = "S2">Set the injection and production flow rates at the target cells:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Cell-centered injection/production flow rates</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">Total      = 150/3600;                                   </span><span class = "S9">% /h &gt;&gt; /s</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw         = zeros(N,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw([1 N])  = [-Total/4 -Total/4];</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw(Nx)     = -Total/4;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw(N-Nx+1) = -Total/4;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">index      = 5101;                                       </span><span class = "S9">% index of center well cell</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Qw(index)  = Total;</span></p></div></div><p class = "S12"><span class = "S2">Next, initialize the fluids properties for the saline water and carbon dioxide phases: </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% fluid properties units: density (Kg/m^3), viscosity (Pa.s), residual saturation (-)</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">water = Fluid(</span><span class = "S13">'saltwater'</span><span class = "S11">,      [1030, 5e-4,   0.3]);   </span><span class = "S9">% resident fluid</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">co2   = Fluid(</span><span class = "S13">'carbon dioxide'</span><span class = "S11">, [950,  7.7e-5, 0.3]);   </span><span class = "S9">% injected fluid</span></p></div></div><p class = "S12"><span class = "S2">We set the initial pressure to 300 bars everywhere in the reservoir and the initial non-wetting phase (CO2) saturation to its residual value:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">P = 300e5.*ones(Nx,Ny,Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">S = co2.Sr.*ones(N,1);</span></p></div></div><p class = "S12"><span class = "S2">Let's add a group of endogeneous (i.e. external) non-wetting CO2 particles which will be injected with the CO2 phase meaning that the injected phase is composed from a mixture of CO2 fluid and particulate suspensions present in the gas stream as a result of the CO2 capture process which cannot be filtered for economic reasons.</span></p><p class = "S3"><span class = "S2">Herein, we assue that salinity induced modibilization of fines is negligeable, hence by leting the colloidal mobilization rate equal to zero the process is effectively excluded as desired. So the processes which are effectively take into account during this simulation are:</span></p><ul class = "S14"><li class = "S15"><span class = "S0">Hydrodynamic mobilization of fine particles for which kinetic rate equals 3.8e-4 m^-1 and the critical fluid velocity amounts to 0.2e-4 m/s</span></li><li class = "S15"><span class = "S0">Deposition into pore bodies with a kinetic rate constant equal to 1.2e-4 m^-1 </span></li><li class = "S15"><span class = "S0">Deposition into pore throats with a kinetic rate constant equal to 6.2e-6 m^-1 </span></li></ul><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">pt1 = Particle(co2, </span><span class = "S16">...</span><span class = "S9">             % fluid phase in which particles flow</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">               [2e-6, </span><span class = "S16">...</span><span class = "S9">           % Mean-size diameter =&gt; for particle population</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                2500, </span><span class = "S16">...</span><span class = "S9">           % Particles density</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                1e-9, </span><span class = "S16">...</span><span class = "S9">           % Diffusion coeff =&gt; important for small size particles</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                0.2e-4, </span><span class = "S16">...</span><span class = "S9">         % Critical fluid velocity</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                0.1, </span><span class = "S16">...</span><span class = "S9">            % Critical fluid salinity</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                3.8e-4, </span><span class = "S16">...</span><span class = "S9">         % Hydrodynamic release rate</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                0, </span><span class = "S16">...</span><span class = "S9">              % Colloidal mobilisation rate</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                1.2e-4, </span><span class = "S16">...</span><span class = "S9">         % Deposition rate in pore surfaces</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                6.2e-6, </span><span class = "S16">...</span><span class = "S9">         % Deposition rate in pore throats</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">                0.6] </span><span class = "S16">...</span><span class = "S9">            % For civan k_phi model</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">                );</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div></div><p class = "S12"><span class = "S2">Initialize different concentrations arrays for the particulate suspensions including the mobile concentration, pore surface deposited concentration, pore throats deposited concentration.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Initialisation of concentration arrays =&gt; particle 1</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">pt1.C     = zeros(N,1);            </span><span class = "S9">% Initial mobile concentration</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">pt1.C_dep = zeros(N,1);            </span><span class = "S9">% Initial surf. deposited conc</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">pt1.C_pt  = zeros(N,1);            </span><span class = "S9">% Initial throats deposited conc</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">pt1.C_tsf = zeros(N,1);            </span><span class = "S9">% Initial mass tranfer &lt;=&gt; phases</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">Cs = zeros(N,1);                   </span><span class = "S9">% Initial salt concentration</span></p></div></div><p class = "S12"><span class = "S2">Make a copy of the initial permeability and porosity fields as they will be updated by the simulator:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% Copy of initial porosity &amp; permeability</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.por0 = Grid.por;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">Grid.K0   = Grid.K;</span></p></div></div><p class = "S12"><span class = "S2">Reset default options of the Newton-Raphson solver to solve the two-phase flow sub-problem. tol is the absolute tolerance for NR convergence, and maxiter is the maximum number of allowed NR iterations befor decision on divergence and local time step subdivision. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">opt.tol     = 1e-6;</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">opt.maxiter = 50;</span></p></div></div><p class = "S12"><span class = "S2">Finally, set the total period of simulation, the number of time steps and the time step used for the sequential simulation algorithm between the two-phase flow sub-system and the particulate transport-kinetics sub-system:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">day = 3600*24;                     </span><span class = "S9">% seconds/day</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">T   = 30*day;                      </span><span class = "S9">% 30 days</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">nt  = 30;                          </span><span class = "S9">% number of time steps</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">dt  = T/nt;                        </span><span class = "S9">% time step</span></p></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S5"><span class = "S2">Sequential simulation loop </span></h2><p class = "S3"><span class = "S2">Now that all the required data and simulation options are set, we're ready to start the simulation of two-phase flow coupled to particle transport. The main simulation loop involved mainly four steps for each global time step: </span></p><ul class = "S14"><li class = "S15"><span class = "S0">First, the two-phase flow (of injected CO2 and resident water) pressure solver is called.The cell-centered pressure is updated after evaluation of each fluid mobility using the non-wetting phase saturation from the previous time step </span></li><li class = "S15"><span class = "S0">Second, the two-phase flow saturation solver is called to update the distribution of the non-wetting phase saturation during this time step. Unlike the pressure solver which is linear the saturation equation is nonlinear. The numerical solution uses an algorithm which repeatedely subdivise the initial time step into internal sub-steps if the most inner Newton-Raphson (NR) iteration loop does not converge during the prescribed number of iterations. This adaptive algorithm is robust and efficient since the Jacobian of the minimized  residuals during the NR iteration is evaluated analytically. </span></li><li class = "S15"><span class = "S0">Once one pressure-saturation are solved, we transport the particles which involves splitting of their convective, kinetic, and diffusive solves.</span></li><li class = "S15"><span class = "S0">Finally, porosity and subsequent permeability change in each grid cell are evaluated and are both updated for the next time step. </span></li></ul><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">tic;              </span><span class = "S9">% start timing the time marching loop</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S16">for </span><span class = "S10">t=1:nt</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   fprintf(</span><span class = "S13">'\n'</span><span class = "S10">); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   fprintf(</span><span class = "S13">'Solving two-phase flow problem. Time = %f days\n'</span><span class = "S10">, t*dt/day);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S9">% call TPFA flow solver adapted to two-phase flow </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   [P,V] = TwoPhasePressure(Grid,S,co2,water,Qw,P,dt); </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S9">% solve for brine saturation </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   S = ImplicitSaturation(Grid,S,co2,water,V,Qw,dt,opt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.sat = reshape(S,Grid.Nx,Grid.Ny,Grid.Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S9">% particles transport solve </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   [m1,~] = RelativePerm(S,co2,water);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   [~,V1] = TwoPhasePressure(Grid,m1.*S,co2,water,Qw);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   pt1 = pt1.transport(Grid,V1,Qw,Cs,dt);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S9">% evaluate porosity and permeability change </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.por0 = reshape(Grid.por0,N,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.por  = reshape(Grid.por,N,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.por  = Grid.por0 - ( (pt1.C_dep + pt1.C_pt) / pt1.density);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   f = abs(ones(N,1) - pt1.alpha_fe.*pt1.C_pt); </span><span class = "S9">% distributed flow efficiency factor</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   k_ratio = EvalPermeabilityCivan(Grid.por0,Grid.por,0,f,3);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.por0 = reshape(Grid.por0,Grid.Nx,Grid.Ny,Grid.Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   Grid.por  = reshape(Grid.por,Grid.Nx,Grid.Ny,Grid.Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   k_ratio = reshape(k_ratio,1,Grid.Nx,Grid.Ny,Grid.Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S16">for </span><span class = "S10">l=1:3</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">       Grid.K(l,:,:,:) = Grid.K0(l,:,:,:) .* k_ratio(1,:,:,:);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">   </span><span class = "S18">end</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">   k_ratio = reshape(k_ratio,N,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S18">end</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 1.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">......<br>Converged in 7 time sub-steps and 221 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.937753e-16.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 2.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">..<br>Converged in 3 time sub-steps and 21 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 3.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">.<br>Converged in 2 time sub-steps and 13 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 4.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 16 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  6.221448e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 5.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 7 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.418084e-16.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 6.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">..<br>Converged in 3 time sub-steps and 19 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 7.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 8 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 8.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">.<br>Converged in 2 time sub-steps and 10 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 9.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 8 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 10.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 8 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 11.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 8 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 12.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 8 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 13.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 9 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 14.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 15 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 15.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 5 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 16.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 5 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 17.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  1.561460e-17.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 18.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 19.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 20.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 21.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 22.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 23.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 24.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 4 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 25.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 26.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsWarningElement" data-width="889" style="width: 889px;"><div class="diagnosticMessage-wrapper diagnosticMessage-warningType"><div class="diagnosticMessage-messagePart">Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND =  7.076031e-19.</div><div class="diagnosticMessage-stackPart"></div></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 27.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 28.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 29.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Solving two-phase flow problem. Time = 30.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Converged in 1 time sub-steps and 3 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Newton-Raphson iteration converged in 1 steps</div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S10">toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="889" style="width: 889px;"><div class="textElement">Elapsed time is 1459.335985 seconds.</div></div></div></div></div><p class = "S12"><span class = "S2">During each time step the output log reports the number of sub-steps and total NR iterations spent to solve the nonlinear saturation equation. Note that each NR iteration corresponds to one call to the inner sparse linear solver whose size equals the number of gid cells. For instance we use the MATLAB backslash </span><span class = "S6">\</span><span class = "S2"> operator (i.e. </span><span class = "S6">mldivide</span><span class = "S2"> direct solver) which is based on </span><a href = 'http://faculty.cse.tamu.edu/davis/suitesparse.html'><span class = "S0">UMFPACK</span></a><span class = "S2"> to solve the inner linear systems. In later versions of this toolkit we plan to introduce iterative Krylov linear solvers with preconditioning which are more efficient for large scale problems. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S5"><span class = "S2">Results visualization and analysis </span></h2><p class = "S3"><span class = "S2"> One important benefit to use a scripting language for numerical simulation tasks is to break the boundaries between model pre-processing, processing, and post-processing. </span></p><p class = "S3"><span class = "S2">In the first stage of this tutorial all needed data was prepared in the same script before going into model execution. This task could be achieved with more efficiency even for reservoirs with much complex geometry and distributed materials and fluids properties. </span></p><p class = "S3"><span class = "S2">Now, we will show how to immediately plot the obtained results at the simulation end without going into storing all computational results into files in open or proprietary formats for post-processing by a third party visualization package. We will plot spatial distributions of the reservoir presuure, P, the carbon dioxide saturation, S, the concentrations of mobile, and pore bodies/throats deposits, and the permeability reduction factor K/K0 (where K0 is the initial reservoir permeability). Herein, we simply use the </span><span class = "S6">contourf</span><span class = "S2"> command to plot filled contour lines of each variable. All subplots are arranged in a 3 x 2 matrix using the </span><span class = "S6">subplot</span><span class = "S2"> command for reporting compactness and each sub-figure has its own title. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot fluid pressure </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,1);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(P,Nx,Ny)'*1e-5,11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'Pressure (bars)'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot CO2 saturation </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,2);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(S,Nx,Ny)',11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'CO_2 saturation'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% mobile co2 particles </span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,3);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(pt1.C,Nx,Ny)',11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'mobile particle concentration'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% deposited co2 particles on pore bodies</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,4);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(pt1.C_dep,Nx,Ny)',11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'pore-body deposited concentration'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% deposited co2 particles on pore throats</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,5);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(pt1.C_pt,Nx,Ny)',11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'pore-throat deposited concentration'</span><span class = "S10">);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10"></span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S9">% plot permeability reduction factor</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">subplot(3,2,6);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S11">         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),</span><span class = "S18">...</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">         reshape(k_ratio,Nx,Ny)',11);</span></p></div><div class = 'inlineWrapper'><p class = "S8 lineNode"><span class = "S10">colormap(jet(16));</span></p></div><div class = 'inlineWrapper outputs'><p class = "S8 lineNode"><span class = "S11">axis </span><span class = "S13">tight equal</span><span class = "S11">; colorbar, title(</span><span class = "S13">'Permeability reduction factor'</span><span class = "S10">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 889px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S12"><span class = "S2">The obtained results reproduce exactly those presented in Figure 6 of Sbai and Azaroual (2011) paper. Mobile particles concentration is most important in a circular zone arround the injection well where CO2 velocity fraction exceeds the prescribed critical velocity leading to hydrodynamic release of particles already deposited in this area. Pore body deposited concentration distribution shape closely follows that of the mobile concentration but with a reversed gardient which is much higher arround the injector. Pore throat deposits and the permeability reduction factor are more restricted to a small radial zone arround the injection well.   </span></p></div></div>
<!-- 
##### SOURCE BEGIN #####
%% Fines migration, hydrodynamic release, and deposition (formation damage) in a five-spot well pattern - case 1
% In this demonstration example we will reproduce the two-dimensional problem 
% for formation damage caused by massive injection of a CO2-particle mixture as 
% discussed by *Sbai and Azaroual (2011) *in their example described in section 
% 5.1. This involves pumpage of a salt aqueous brine from a 3000m deep reservoir 
% through four square corner's wells with an equal injection rate. CO2 phase is 
% injected from a well at the center of the domain with a total injection flow 
% rate of 1.15 x 10^6 tons/year. 
%% Reservoir, fluids,and particulate suspensions properties setup 
% Let's first specify the simple reservoir geometry which is a one layer of 
% 500m x 500m x 5m along each space dimension. Additionally, we will explicitly 
% store some required variables in the '|Grid|' structure as reequired for next 
% computational tasks. These are: |Nx|, |Ny|, |Nz| which corresponds to the number 
% of cells along each space dismension, respectively. N: the total number of grid 
% cells, |hx|, |hy|, |hz|, which are the uniform spacing along each space dimension, 
% |V|: the volume of grid cells, |K|: the initial distribution of the permeability, 
% |compr| the total compressibility factor, |por|: the initial distribution of 
% the porosity. 

% Domain size along X, Y & Z directions
dom = RecDomain([500,500,5]);
Dx = dom.Lx; Dy = dom.Ly; Dz = dom.Lz;
Grid.Nx = 101;  Grid.Ny = 101;   Grid.Nz = 1; 
Nx = Grid.Nx;   Ny = Grid.Ny;    Nz = Grid.Nz; 
Grid.hx = (Dx/Nx);
Grid.hy = (Dy/Ny);
Grid.hz = (Dz/Nz);
Grid.N = Grid.Nx*Grid.Ny*Grid.Nz;
N = Grid.N;
Grid.V = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);                        % grid cells volumes
Grid.K = 0.85e-12.*ones(3,Nx,Ny,Nz);                      
Grid.compr = 4.4e-10.*ones(Grid.Nx,Grid.Ny,Grid.Nz);     % compressibility
Grid.por   = 0.3.*ones(Grid.Nx,Grid.Ny,Grid.Nz);         % porosity

%% 
% Set the injection and production flow rates at the target cells:

% Cell-centered injection/production flow rates
Total      = 150/3600;                                   % /h >> /s
Qw         = zeros(N,1);
Qw([1 N])  = [-Total/4 -Total/4];
Qw(Nx)     = -Total/4;
Qw(N-Nx+1) = -Total/4;
index      = 5101;                                       % index of center well cell
Qw(index)  = Total;
%% 
% Next, initialize the fluids properties for the saline water and carbon 
% dioxide phases: 

% fluid properties units: density (Kg/m^3), viscosity (Pa.s), residual saturation (-)
water = Fluid('saltwater',      [1030, 5e-4,   0.3]);   % resident fluid
co2   = Fluid('carbon dioxide', [950,  7.7e-5, 0.3]);   % injected fluid
%% 
% We set the initial pressure to 300 bars everywhere in the reservoir and 
% the initial non-wetting phase (CO2) saturation to its residual value:

P = 300e5.*ones(Nx,Ny,Nz);
S = co2.Sr.*ones(N,1);
%% 
% Let's add a group of endogeneous (i.e. external) non-wetting CO2 particles 
% which will be injected with the CO2 phase meaning that the injected phase is 
% composed from a mixture of CO2 fluid and particulate suspensions present in 
% the gas stream as a result of the CO2 capture process which cannot be filtered 
% for economic reasons.
% 
% Herein, we assue that salinity induced modibilization of fines is negligeable, 
% hence by leting the colloidal mobilization rate equal to zero the process is 
% effectively excluded as desired. So the processes which are effectively take 
% into account during this simulation are:
% 
% * Hydrodynamic mobilization of fine particles for which kinetic rate equals 
% 3.8e-4 m^-1 and the critical fluid velocity amounts to 0.2e-4 m/s
% * Deposition into pore bodies with a kinetic rate constant equal to 1.2e-4 
% m^-1 
% * Deposition into pore throats with a kinetic rate constant equal to 6.2e-6 
% m^-1 

pt1 = Particle(co2, ...             % fluid phase in which particles flow
               [2e-6, ...           % Mean-size diameter => for particle population
                2500, ...           % Particles density
                1e-9, ...           % Diffusion coeff => important for small size particles
                0.2e-4, ...         % Critical fluid velocity
                0.1, ...            % Critical fluid salinity
                3.8e-4, ...         % Hydrodynamic release rate
                0, ...              % Colloidal mobilisation rate
                1.2e-4, ...         % Deposition rate in pore surfaces
                6.2e-6, ...         % Deposition rate in pore throats
                0.6] ...            % For civan k_phi model
                );

%% 
% Initialize different concentrations arrays for the particulate suspensions 
% including the mobile concentration, pore surface deposited concentration, pore 
% throats deposited concentration.

% Initialisation of concentration arrays => particle 1
pt1.C     = zeros(N,1);            % Initial mobile concentration
pt1.C_dep = zeros(N,1);            % Initial surf. deposited conc
pt1.C_pt  = zeros(N,1);            % Initial throats deposited conc
pt1.C_tsf = zeros(N,1);            % Initial mass tranfer <=> phases
Cs = zeros(N,1);                   % Initial salt concentration
%% 
% Make a copy of the initial permeability and porosity fields as they will 
% be updated by the simulator:

% Copy of initial porosity & permeability
Grid.por0 = Grid.por;
Grid.K0   = Grid.K;
%% 
% Reset default options of the Newton-Raphson solver to solve the two-phase 
% flow sub-problem. tol is the absolute tolerance for NR convergence, and maxiter 
% is the maximum number of allowed NR iterations befor decision on divergence 
% and local time step subdivision. 

opt.tol     = 1e-6;
opt.maxiter = 50;
%% 
% Finally, set the total period of simulation, the number of time steps 
% and the time step used for the sequential simulation algorithm between the two-phase 
% flow sub-system and the particulate transport-kinetics sub-system:

day = 3600*24;                     % seconds/day
T   = 30*day;                      % 30 days
nt  = 30;                          % number of time steps
dt  = T/nt;                        % time step
%% Sequential simulation loop 
% Now that all the required data and simulation options are set, we're ready 
% to start the simulation of two-phase flow coupled to particle transport. The 
% main simulation loop involved mainly four steps for each global time step: 
% 
% * First, the two-phase flow (of injected CO2 and resident water) pressure 
% solver is called.The cell-centered pressure is updated after evaluation of each 
% fluid mobility using the non-wetting phase saturation from the previous time 
% step 
% * Second, the two-phase flow saturation solver is called to update the distribution 
% of the non-wetting phase saturation during this time step. Unlike the pressure 
% solver which is linear the saturation equation is nonlinear. The numerical solution 
% uses an algorithm which repeatedely subdivise the initial time step into internal 
% sub-steps if the most inner Newton-Raphson (NR) iteration loop does not converge 
% during the prescribed number of iterations. This adaptive algorithm is robust 
% and efficient since the Jacobian of the minimized  residuals during the NR iteration 
% is evaluated analytically. 
% * Once one pressure-saturation are solved, we transport the particles which 
% involves splitting of their convective, kinetic, and diffusive solves.
% * Finally, porosity and subsequent permeability change in each grid cell are 
% evaluated and are both updated for the next time step. 

tic;              % start timing the time marching loop
for t=1:nt
   fprintf('\n'); 
   fprintf('Solving two-phase flow problem. Time = %f days\n', t*dt/day);
   
   % call TPFA flow solver adapted to two-phase flow 
   [P,V] = TwoPhasePressure(Grid,S,co2,water,Qw,P,dt); 
   
   % solve for brine saturation 
   S = ImplicitSaturation(Grid,S,co2,water,V,Qw,dt,opt);
   Grid.sat = reshape(S,Grid.Nx,Grid.Ny,Grid.Nz);
   
   % particles transport solve 
   [m1,~] = RelativePerm(S,co2,water);
   [~,V1] = TwoPhasePressure(Grid,m1.*S,co2,water,Qw);
   pt1 = pt1.transport(Grid,V1,Qw,Cs,dt);
   
   % evaluate porosity and permeability change 
   Grid.por0 = reshape(Grid.por0,N,1);
   Grid.por  = reshape(Grid.por,N,1);
   Grid.por  = Grid.por0 - ( (pt1.C_dep + pt1.C_pt) / pt1.density);
   f = abs(ones(N,1) - pt1.alpha_fe.*pt1.C_pt); % distributed flow efficiency factor
   k_ratio = EvalPermeabilityCivan(Grid.por0,Grid.por,0,f,3);
   Grid.por0 = reshape(Grid.por0,Grid.Nx,Grid.Ny,Grid.Nz);
   Grid.por  = reshape(Grid.por,Grid.Nx,Grid.Ny,Grid.Nz);
   
   k_ratio = reshape(k_ratio,1,Grid.Nx,Grid.Ny,Grid.Nz);
   for l=1:3
       Grid.K(l,:,:,:) = Grid.K0(l,:,:,:) .* k_ratio(1,:,:,:);
   end
   k_ratio = reshape(k_ratio,N,1);

end
toc;
%% 
% During each time step the output log reports the number of sub-steps and 
% total NR iterations spent to solve the nonlinear saturation equation. Note that 
% each NR iteration corresponds to one call to the inner sparse linear solver 
% whose size equals the number of gid cells. For instance we use the MATLAB backslash 
% |\| operator (i.e. |mldivide| direct solver) which is based on <http://faculty.cse.tamu.edu/davis/suitesparse.html 
% UMFPACK> to solve the inner linear systems. In later versions of this toolkit 
% we plan to introduce iterative Krylov linear solvers with preconditioning which 
% are more efficient for large scale problems. 
%% Results visualization and analysis 
%  One important benefit to use a scripting language for numerical simulation 
% tasks is to break the boundaries between model pre-processing, processing, and 
% post-processing. 
% 
% In the first stage of this tutorial all needed data was prepared in the 
% same script before going into model execution. This task could be achieved with 
% more efficiency even for reservoirs with much complex geometry and distributed 
% materials and fluids properties. 
% 
% Now, we will show how to immediately plot the obtained results at the simulation 
% end without going into storing all computational results into files in open 
% or proprietary formats for post-processing by a third party visualization package. 
% We will plot spatial distributions of the reservoir presuure, P, the carbon 
% dioxide saturation, S, the concentrations of mobile, and pore bodies/throats 
% deposits, and the permeability reduction factor K/K0 (where K0 is the initial 
% reservoir permeability). Herein, we simply use the |contourf| command to plot 
% filled contour lines of each variable. All subplots are arranged in a 3 x 2 
% matrix using the |subplot| command for reporting compactness and each sub-figure 
% has its own title. 

figure; 

% plot fluid pressure 
subplot(3,2,1);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(P,Nx,Ny)'*1e-5,11);
colormap(jet(16));
axis tight equal; colorbar, title('Pressure (bars)');

% plot CO2 saturation 
subplot(3,2,2);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(S,Nx,Ny)',11);
colormap(jet(16));
axis tight equal; colorbar, title('CO_2 saturation');

% mobile co2 particles 
subplot(3,2,3);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(pt1.C,Nx,Ny)',11);
colormap(jet(16));
axis tight equal; colorbar, title('mobile particle concentration');

% deposited co2 particles on pore bodies
subplot(3,2,4);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(pt1.C_dep,Nx,Ny)',11);
colormap(jet(16));
axis tight equal; colorbar, title('pore-body deposited concentration');

% deposited co2 particles on pore throats
subplot(3,2,5);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(pt1.C_pt,Nx,Ny)',11);
colormap(jet(16));
axis tight equal; colorbar, title('pore-throat deposited concentration');

% plot permeability reduction factor
subplot(3,2,6);
contourf(linspace((Dx/Nx)/2,Dx-(Dx/Nx)/2,Nx),...
         linspace((Dy/Ny)/2,Dy-(Dy/Ny)/2,Ny),...
         reshape(k_ratio,Nx,Ny)',11);
colormap(jet(16));
axis tight equal; colorbar, title('Permeability reduction factor');
%% 
% The obtained results reproduce exactly those presented in Figure 6 of 
% Sbai and Azaroual (2011) paper. Mobile particles concentration is most important 
% in a circular zone arround the injection well where CO2 velocity fraction exceeds 
% the prescribed critical velocity leading to hydrodynamic release of particles 
% already deposited in this area. Pore body deposited concentration distribution 
% shape closely follows that of the mobile concentration but with a reversed gardient 
% which is much higher arround the injector. Pore throat deposits and the permeability 
% reduction factor are more restricted to a small radial zone arround the injection 
% well.
##### SOURCE END #####
--></body></html>